{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e4c28a13-49eb-4b04-8e96-b62f001e8406",
   "metadata": {},
   "source": [
    "# 最大似然估计\n",
    "\n",
    "原文: https://blog.csdn.net/lly1122334/article/details/105046070"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "399f2dbc-7b1d-4091-9ef2-387f1de03d14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAARU0lEQVR4nO3df4zkdX3H8edLwN+mQDmv593Ro4oaNAXMFrG2DUpVRONhYgmmxVNpzjRotSFVsEm1aWmurUpptDSnUI6WggSxXCy1UqQ1JhU8EPl1Wq8KcteDOwSR1pT28N0/5nsyLru3P2ZnZ/azz0eyme98vt+Zee/d7mvf85nPfCdVhSSpLU8ZdQGSpIVnuEtSgwx3SWqQ4S5JDTLcJalBB4+6AIAjjjii1q1bN+oyJGlJueWWWx6sqhVT7RuLcF+3bh3btm0bdRmStKQkuXe6fU7LSFKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg8biHarSUrfu3H/48fY9m94wwkqkHjt3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KAZwz3J05PcnOTrSe5K8gfd+FFJbkqyI8mnkzy1G39ad31Ht3/dkL8HSdIks+ncHwNeXVXHAscBpyQ5EfgT4IKqegHwMHBWd/xZwMPd+AXdcZKkRTRjuFfPf3VXD+m+Cng1cHU3vgU4rdte312n239ykixUwZKkmc3qHapJDgJuAV4AfAL4D+D7VbWvO2QnsLrbXg3cB1BV+5I8Avw08OCk+9wIbAQ48sgjB/supEUyyDtRfRerFtOsXlCtqser6jhgDXAC8OJBH7iqNlfVRFVNrFgx5Yd3S5LmaU7nlqmq7ye5EXgFcGiSg7vufQ2wqztsF7AW2JnkYOCngO8tYM3SkmcXr2GbzWqZFUkO7bafAbwG2A7cCLylO2wDcG23vbW7Trf/i1VVC1izJGkGs+ncVwFbunn3pwBXVdXnktwNXJnkj4CvARd3x18M/E2SHcBDwBlDqFuSdAAzhntV3Q4cP8X4t+nNv08e/x/g1xakOmmM9U+tDHLMdMc7XaNB+A5VSWqQ4S5JDTLcJalBhrskNchwl6QG+QHZUseVKmqJnbskNchwl6QGGe6S1CDDXZIaZLhLUoNcLSPNYK7nh5HGgZ27JDXIcJekBhnuktQg59ylKQwyzz6MOXrfPau5MtylEfMDPTQMTstIUoMMd0lqkOEuSQ1yzl3Nc45ay5GduyQ1yHCXpAY5LSONKc9po0HM2LknWZvkxiR3J7kryXu78Q8n2ZXktu7r1L7bnJdkR5JvJnndML8BSdKTzaZz3wecU1W3JnkOcEuS67t9F1TVR/oPTnIMcAbwEuB5wD8neWFVPb6QhUuDsjNWy2bs3Ktqd1Xd2m0/CmwHVh/gJuuBK6vqsar6DrADOGEhipUkzc6cXlBNsg44HripG3p3ktuTXJLksG5sNXBf3812MsUfgyQbk2xLsm3v3r1zr1ySNK1Zh3uSZwOfAd5XVT8ALgKeDxwH7AY+OpcHrqrNVTVRVRMrVqyYy00lSTOYVbgnOYResF9eVdcAVNUDVfV4Vf0I+CRPTL3sAtb23XxNNyZJWiSzWS0T4GJge1V9rG98Vd9hbwbu7La3AmckeVqSo4CjgZsXrmRJ0kxms1rmlcCZwB1JbuvGPgi8NclxQAH3AO8CqKq7klwF3E1vpc3ZrpSRpMU1Y7hX1ZeBTLHrugPc5nzg/AHqkiQNwNMPSFKDDHdJapDhLkkN8sRh0hLmueo1HcNdTZruvDGeT0bLhdMyktQgw12SGuS0jJrhlIv0BDt3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoN8E5PUCE8ipn527pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNWjGde5J1gKXASuBAjZX1YVJDgc+DawD7gFOr6qHkwS4EDgV+CHw9qq6dTjlazlyPbc0s9l07vuAc6rqGOBE4OwkxwDnAjdU1dHADd11gNcDR3dfG4GLFrxqSdIBzRjuVbV7f+ddVY8C24HVwHpgS3fYFuC0bns9cFn1fAU4NMmqhS5ckjS9Oc25J1kHHA/cBKysqt3drvvpTdtAL/jv67vZzm5s8n1tTLItyba9e/fOtW5J0gHMOtyTPBv4DPC+qvpB/76qKnrz8bNWVZuraqKqJlasWDGXm0qSZjCrE4clOYResF9eVdd0ww8kWVVVu7tplz3d+C5gbd/N13RjkhaJLzprxs69W/1yMbC9qj7Wt2srsKHb3gBc2zf+tvScCDzSN30jSVoEs+ncXwmcCdyR5LZu7IPAJuCqJGcB9wKnd/uuo7cMcge9pZDvWMiCJUkzmzHcq+rLQKbZffIUxxdw9oB1SVogTtEsT35Yh7SMGPTLh6cfkKQGGe6S1CDDXZIaZLhLUoN8QVVapnxxtW127pLUIDt3LWn93aekJ9i5S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg3yHqsaK5zuRFoaduyQ1yHCXpAY5LaOx5RSNNH927pLUIMNdkho0Y7gnuSTJniR39o19OMmuJLd1X6f27TsvyY4k30zyumEVLkma3mw690uBU6YYv6Cqjuu+rgNIcgxwBvCS7jZ/meSghSpWkjQ7M4Z7VX0JeGiW97ceuLKqHquq7wA7gBMGqE+SNA+DzLm/O8nt3bTNYd3YauC+vmN2dmNPkmRjkm1Jtu3du3eAMiRJk813KeRFwB8C1V1+FHjnXO6gqjYDmwEmJiZqnnVomfCzUkfHJalL07zCvaoe2L+d5JPA57qru4C1fYeu6cYkjTH/eLZnXtMySVb1XX0zsH8lzVbgjCRPS3IUcDRw82AlSpLmasbOPckVwEnAEUl2Ah8CTkpyHL1pmXuAdwFU1V1JrgLuBvYBZ1fV40OpXJI0rRnDvareOsXwxQc4/nzg/EGKkiQNxnPLaOSc75UWnqcfkKQGGe6S1CDDXZIaZLhLUoMMd0lqkKtlNBKukJGGy85dkhpkuEtSgwx3SWqQc+6SZs3T/y4ddu6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNmjHck1ySZE+SO/vGDk9yfZJvdZeHdeNJ8hdJdiS5PcnLhlm8JGlqsznl76XAx4HL+sbOBW6oqk1Jzu2ufwB4PXB09/Vy4KLuUvKj9Rrj6X/H24yde1V9CXho0vB6YEu3vQU4rW/8sur5CnBoklULVKskaZbm+2EdK6tqd7d9P7Cy214N3Nd33M5ubDeTJNkIbAQ48sgj51mGpHFgFz9+Bn5BtaoKqHncbnNVTVTVxIoVKwYtQ5LUZ77h/sD+6Zbuck83vgtY23fcmm5MkrSI5hvuW4EN3fYG4Nq+8bd1q2ZOBB7pm76RJC2SGefck1wBnAQckWQn8CFgE3BVkrOAe4HTu8OvA04FdgA/BN4xhJolSTOYMdyr6q3T7Dp5imMLOHvQoiRJg5nvahlpWq6ckEbPcJc0NP6hHx3PLSNJDbJz11B5ygFpNOzcJalBhrskNchwl6QGGe6S1CDDXZIa5GoZSQtqNiukXP8+fHbuktQgw12SGmS4S1KDnHPXgvCdqNJ4sXOXpAYZ7pLUIMNdkhrknLukReHrMovLcNe8+csqjS+nZSSpQYa7JDXIcJekBjnnrjlxnl1aGuzcJalBA3XuSe4BHgUeB/ZV1USSw4FPA+uAe4DTq+rhwcqUtNx4WuDBLMS0zKuq6sG+6+cCN1TVpiTndtc/sACPoxFxKkZaeoYx574eOKnb3gL8C4a7pGnYoQ/HoHPuBXwhyS1JNnZjK6tqd7d9P7Byqhsm2ZhkW5Jte/fuHbAMSVK/QTv3X6qqXUmeC1yf5Bv9O6uqktRUN6yqzcBmgImJiSmPkbS8OAW4cAbq3KtqV3e5B/gscALwQJJVAN3lnkGLlCTNzbzDPcmzkjxn/zbwWuBOYCuwoTtsA3DtoEVKkuZmkGmZlcBnk+y/n7+rqs8n+SpwVZKzgHuB0wcvU4vBF7akdsw73Kvq28CxU4x/Dzh5kKIkSYPxHaqS1CDPLaMpuWpBWtrs3CWpQYa7JDXIcJekBhnuktQgw12SGuRqGUljzzfYzZ2duyQ1yM59GbIL0lLmz+/sGO6Nmu0vgG9WktrktIwkNcjOfRmwO9dy49SNnbskNcnOXVITfIb6k+zcJalBdu6Sliy79enZuUtSgwx3SWqQ0zKSmrZcl0Ua7kvQdPOMy+kHV5qP5RT0hrskzcLkpmrc/zgY7kvEbFYFuHJAmr3WnwEb7mPMsJY0X0ML9ySnABcCBwGfqqpNw3qscTCbubzpjllO84DSuGulqRpKuCc5CPgE8BpgJ/DVJFur6u5hPN5SNN0PUCs/WFLrZjOtM8rGbVid+wnAjqr6NkCSK4H1wIKH+0L+4831vpwHlzTZXBu3YYV+qmrh7zR5C3BKVf1md/1M4OVV9e6+YzYCG7urLwK+uUAPfwTw4ALd17CMe43jXh+Mf43jXh+Mf43jXh+MvsafraoVU+0Y2QuqVbUZ2LzQ95tkW1VNLPT9LqRxr3Hc64Pxr3Hc64Pxr3Hc64PxrnFYpx/YBaztu76mG5MkLYJhhftXgaOTHJXkqcAZwNYhPZYkaZKhTMtU1b4k7wb+id5SyEuq6q5hPNYUFnyqZwjGvcZxrw/Gv8Zxrw/Gv8Zxrw/GuMahvKAqSRotT/krSQ0y3CWpQc2Ee5K1SW5McneSu5K8d9Q1TSXJQUm+luRzo65lKkkOTXJ1km8k2Z7kFaOuqV+S3+n+f+9MckWSp49BTZck2ZPkzr6xw5Ncn+Rb3eVhY1jjn3X/z7cn+WySQ8epvr595ySpJEeMora+OqasMcl7un/Hu5L86ajqm6yZcAf2AedU1THAicDZSY4ZcU1TeS+wfdRFHMCFwOer6sXAsYxRrUlWA78NTFTVS+m9WH/GaKsC4FLglElj5wI3VNXRwA3d9VG6lCfXeD3w0qr6eeDfgfMWu6g+l/Lk+kiyFngt8N3FLmgKlzKpxiSvovfu+2Or6iXAR0ZQ15SaCfeq2l1Vt3bbj9ILpdWjreonJVkDvAH41KhrmUqSnwJ+BbgYoKr+t6q+P9Kinuxg4BlJDgaeCfzniOuhqr4EPDRpeD2wpdveApy2mDVNNlWNVfWFqtrXXf0KvfejjMQ0/4YAFwDvB0a+8mOaGn8L2FRVj3XH7Fn0wqbRTLj3S7IOOB64acSlTPbn9H5QfzTiOqZzFLAX+Otu6uhTSZ416qL2q6pd9Dqj7wK7gUeq6gujrWpaK6tqd7d9P7BylMXMwjuBfxx1Ef2SrAd2VdXXR13LAbwQ+OUkNyX51yS/MOqC9msu3JM8G/gM8L6q+sGo69kvyRuBPVV1y6hrOYCDgZcBF1XV8cB/M/rphB/r5q3X0/sj9DzgWUl+Y7RVzax6641H3nlOJ8nv0ZvWvHzUteyX5JnAB4HfH3UtMzgYOJzeVPDvAlclyWhL6mkq3JMcQi/YL6+qa0ZdzySvBN6U5B7gSuDVSf52tCU9yU5gZ1Xtf8ZzNb2wHxe/CnynqvZW1f8B1wC/OOKapvNAklUA3eXYPF3vl+TtwBuBX6/xetPL8+n9Ef969zuzBrg1yc+MtKon2wlcUz0303tWPtIXfvdrJty7v5YXA9ur6mOjrmeyqjqvqtZU1Tp6LwJ+sarGquusqvuB+5K8qBs6mSGcpnkA3wVOTPLM7v/7ZMboBd9JtgIbuu0NwLUjrGVK3QfqvB94U1X9cNT19KuqO6rquVW1rvud2Qm8rPsZHSd/D7wKIMkLgacyJmeybCbc6XXGZ9LriG/rvk4ddVFL0HuAy5PcDhwH/PFoy3lC94ziauBW4A56P78jf/t3kiuAfwNelGRnkrOATcBrknyL3jOOkX4S2TQ1fhx4DnB99/vyV2NW31iZpsZLgJ/rlkdeCWwYl2dAnn5AkhrUUucuSeoY7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalB/w8/S/LJYDw6wAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10.013448891998049, 1.9957983058594568)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "u = 10 # 数学期望\n",
    "o = 2 #方差\n",
    "x = u + o * np.random.randn(10000) # 正态分布\n",
    "plt.hist(x, bins=100)\n",
    "plt.show()\n",
    "print(norm.fit(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5360864-7b47-4d07-91de-723eb549a02f",
   "metadata": {},
   "source": [
    "## 分析\n",
    "原文: https://cloud.tencent.com/developer/article/1189259\n",
    "\n",
    "极大似然估计（Maximum likelihood estimation, 简称MLE）是很常用的参数估计方法，极大似然原理的直观想法是，一个随机试验如有若干个可能的结果A，B，C，... ，若在一次试验中，结果A出现了，那么可以认为实验条件对A的出现有利，也即出现的概率P(A)较大。也就是说，如果已知某个随机样本满足某种概率分布，但是其中具体的参数不清楚，参数估计就是通过若干次试验，观察其结果，利用结果推出参数的大概值。极大似然估计是建立在这样的思想上：已知某个参数能使这个样本出现的概率最大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值（请参见“百度百科”）。\n",
    "\n",
    "本文以一个简单的离散型分布的例子，模拟投掷硬币估计头像(head)向上的概率。投掷硬币落到地面后，不是head向上就是tail朝上，这是一个典型的伯努利实验，形成一个伯努利分布，有着如下的离散概率分布函数：\n",
    "\n",
    "$$ \\phi(x)=p^x(1-p)^{(1-x)} $$\n",
    "\n",
    "其中，x等于1或者0，即结果，这里用1表示head、0表示tail。\n",
    "\n",
    "对于n次独立的投掷，很容易写出其似然函数：\n",
    "\n",
    "$$ \\mathcal{L}(p ; \\mathbf{x})=\\prod_{i=1}^{n} p^{x_{i}}(1-p)^{1-x_{i}} $$\n",
    "\n",
    "现在想用极大似然估计的方法把p估计出来。就是使得上面这个似然函数取极大值的情况下的p的取值，就是要估计的参数。\n",
    "\n",
    "首先用Python把投掷硬币模拟出来："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "0178d1a5-276e-4c71-869e-295529282976",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 1 0 0 0 0 0 1 1 1 0 0 1 0 0 1 0 1 1 0 0 1 1 1 1 1 0 0 1 0]\n",
      "p**51*(1 - p)**49\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import bernoulli\n",
    "p_1 = 1.0/2           # 假设这是我们需要估计的\n",
    "fp = bernoulli(p_1)   # 产生伯努力随机变量\n",
    "xs = fp.rvs(100)      # 产生100个样本\n",
    "print(xs[:30])        # 看看前面30个\n",
    "# 通过此模拟，使用sympy库把似然函数写出来：\n",
    "import sympy\n",
    "import numpy as np\n",
    "x, p, z = sympy.symbols('x p z', positive=True)\n",
    "phi = p**x*(1-p)**(1-x)                          # 分布函数\n",
    "L = np.prod([phi.subs(x, i) for i in xs])        # 似然函数\n",
    "print(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05832dea-f348-4b1d-9f9d-d151dd441306",
   "metadata": {},
   "source": [
    "从上面的结论可以看出，作100次伯努利实验，出现positive、1及head的数目是53个，相应的0也就是tail的数目是47个，比较接近我们设的初始值0.5即1.0/2（注意：现在我们假设p是未知的，要去估计它，看它经过Python的极大似然估计是不是0.5！）。\n",
    "\n",
    "下面，我们使用Python求解这个似然函数取极大值时的p值："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "2dbc9c7f-7d78-41ac-bd8d-a348cc17f254",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[51/100]\n"
     ]
    }
   ],
   "source": [
    "logL = sympy.expand_log(sympy.log(L))\n",
    "sol=sympy.solve(sympy.diff(logL,p),p)\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fffafd27-a216-4258-ae42-68f4c8774658",
   "metadata": {},
   "source": [
    "结果没有什么悬念，53/100的值很接近0.5！\n",
    "\n",
    "取对数后，上面Python的算法最后实际上是求解下式为0的p值：\n",
    "\n",
    "$$ \\log (\\mathcal{L}(p \\mid \\mathbf{x}))=\\log (p) \\sum_{i=1}^{n} x_{i}+\\log (1-p)\\left(n-\\sum_{i=1}^{n} x_{i}\\right) $$\n",
    "\n",
    "上式留给网友自行推导，很多资料都可找到该式。这个式子，是著名的Logistic回归参数估计的极大似然估计算法的基础。\n",
    "\n",
    "进一步，为了更加直观的理解投掷硬币的伯努利实验，我们给出以均值（均值为100*0.5=50）为中心对称的加总离散概率（概率质量函数（probability mass function），Python里面使用pmf函数计算）："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "1ba1daf2-ea3c-4a47-9289-3e9a83037cfb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<scipy.stats._distn_infrastructure.rv_frozen object at 0x7f22b5096670>\n",
      "nan\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import binom\n",
    "b=binom(100, 5)\n",
    "print(b)\n",
    "# 以均值为中心对称的加总概率\n",
    "g = lambda i:b.pmf(np.arange(-i, i) + 50).sum()\n",
    "print(g(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc59d452-7e46-4972-bb70-46e93f687129",
   "metadata": {},
   "source": [
    "本文针对简单的离散概率质量函数的分布使用Python进行了极大似然估计，同时该方法可以应用于连续分布的情形，只要通过其概率密度函数得出其似然函数即可。希望网友把本文的代码实践一遍，也可以和R语言、SAS等软件得到的结论相比较，从而得到更好的极大似然估计的实现方法。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
